#Replication Paper
#Gov 2001
#Anna Hopper & Sam Imlay
#April 27, 2013

setwd("~/Dropbox/ReplicationPaper")

#packages needed
library(survival)
library(vcd)

##########-----------------
#MULTIPLE IMPUTATION

newdata2<-read.csv("~/Dropbox/A-K/Gov2000/imputable.csv")
head(newdata2)
#make partyfeeling numeric 
newdata2$partyfeeling1<-as.numeric(newdata2$partyfeeling1)
newdata2$partyfeeling2<-as.numeric(newdata2$partyfeeling2)
newdata2$partyfeeling3<-as.numeric(newdata2$partyfeeling3)

#redefine votechoice so it can be used for imputation 
newdata2$votechoice<-0
for (i in 1:887){
ifelse(newdata2$votechoice1[i]==1, newdata2$votechoice[i]<-1, newdata2$votechoice[i])
ifelse(newdata2$votechoice2[i]==1, newdata2$votechoice[i]<-2, newdata2$votechoice[i])
ifelse(newdata2$votechoice3[i]==1, newdata2$votechoice[i]<-3, newdata2$votechoice[i])
}

head(newdata2$votechoice)
library(Amelia)
a.out<-amelia(x=newdata2, m=5, noms=c("outsidecon1", "outsidecon2","outsidecon3"), idvars=(c("aconsnam", "votechoice1", "votechoice2", "votechoice3","vote", "newincumb2", "newincumb3", "newincumb1", "labour1", "conservatives1", "liberal_democrats1", "labour2", "conservatives2", "liberal_democrats2", "labour3", "conservatives3", "liberal_democrats3", "touse")))
a.out
save(a.out, file = "imputations.RData")
write.amelia(obj=a.out, file.stem = "imputations")

##-------------
#Since we aren't using Zelig for simulation, we have to combine imputed datasets by simulated 200 betas with each

#First dataset
Idata1<-read.csv("~/Dropbox/a-k/Gov2000/Idata1.csv")
clogitI1<-clogit(votechoice~conservatives+liberal_democrats+partyFeeling+drdistp+inccand+absdepdist+outsideCon+strata(id)+cluster(aconsnam), method="approximate", data=Idata1)
clogitI1
Idata1$partyfeeling<-Idata1$partyFeeling
Idata1$outsidecon<-Idata1$outsideCon
Idata1$partyFeeling<-NULL
Idata1$outsideCon<-NULL
dim(Idata1)

meanI1<-coef(clogitI1)
vcvI1<-vcov(clogitI1)

#randomly draw betas 
library(mvtnorm)
set.seed(02138)
betasI1<-rmvnorm(200, meanI1, vcvI1) 



#Second dataset
Idata2<-read.csv("~/Dropbox/a-k/Gov2000/Idata2.csv")

clogitI2<-clogit(votechoice~conservatives+liberal_democrats+partyfeeling+drdistp+inccand+absdepdist+outsidecon+strata(id)+cluster(aconsnam), method="approximate", data=Idata2)
clogitI2
dim(Idata2)


meanI2<-coef(clogitI2)
vcvI2<-vcov(clogitI2)

#randomly draw betas 
library(mvtnorm)
set.seed(02138)
betasI2<-rmvnorm(200, meanI2, vcvI2) 

##Third dataset
Idata3<-read.csv("~/Dropbox/a-k/Gov2000/Idata3.csv")
clogitI3<-clogit(votechoice~conservatives+liberal_democrats+partyfeeling+drdistp+inccand+absdepdist+outsidecon+strata(id)+cluster(aconsnam), method="approximate", data=Idata3)
clogitI3

meanI3<-coef(clogitI3)
vcvI3<-vcov(clogitI3)

#randomly draw betas 
library(mvtnorm)
set.seed(02138)
betasI3<-rmvnorm(200, meanI3, vcvI3) 

##Fourth dataset
Idata4<-read.csv("~/Dropbox/a-k/Gov2000/Idata4.csv")
clogitI4<-clogit(votechoice~conservatives+liberal_democrats+partyfeeling+drdistp+inccand+absdepdist+outsidecon+strata(id)+cluster(aconsnam), method="approximate", data=Idata4)
clogitI4

meanI4<-coef(clogitI4)
vcvI4<-vcov(clogitI4)

#randomly draw betas 
library(mvtnorm)
set.seed(02138)
betasI4<-rmvnorm(200, meanI4, vcvI4) 

##Fifth dataset 
Idata5<-read.csv("~/Dropbox/a-k/Gov2000/Idata5.csv")
clogitI5<-clogit(votechoice~conservatives+liberal_democrats+partyfeeling+drdistp+inccand+absdepdist+outsidecon+strata(id)+cluster(aconsnam), method="approximate", data=Idata5)
clogitI5

meanI5<-coef(clogitI5)
vcvI5<-vcov(clogitI5)

#randomly draw betas 
library(mvtnorm)
set.seed(02138)
betasI5<-rmvnorm(200, meanI5, vcvI5) 

head(betasI5)

MIbetas<-matrix(nrow=1000, ncol=3)
MIbetas<-rbind(betasI1, betasI2, betasI3, betasI4, betasI5)

IMPUTE<-rbind(Idata1, Idata2, Idata3, Idata4, Idata5)
IMPUTE<-as.data.frame(IMPUTE)

####RESULTS for IMPUTATION
##Subset by Row
LabRows<-subset(IMPUTE, IMPUTE$partyvote==1)

#Conservative matrix: 
ConsRows<-subset(IMPUTE, IMPUTE$partyvote==2)

#Liberal Democrat matrix: 
LDRows<-subset(IMPUTE, IMPUTE$partyvote==3)


#simulate distances 
dists<-seq(0, 120, by=5)

p.ests.cons<-matrix(data=NA, ncol=length(dists), nrow=1000)

for (m in 1:length(dists)){

#empty vectors for pi values to go in 
piL<-NULL
piC<-NULL
piLD<-NULL

# Insert real mean distances for each party
		j<-mean(IMPUTE$drdistp)
		k<-dists[m]
		l<-mean(IMPUTE$drdistp)
	# Establish X values
xL<-c(mean(LabRows$conservatives), mean(LabRows$liberal_democrats), mean(LabRows$partyfeeling), j, mean(LabRows$inccand), median(LabRows$absdepdist), mean(LabRows$outsidecon))

xC<-c(mean(ConsRows$conservatives), mean(ConsRows$liberal_democrats), mean(ConsRows$partyfeeling), k, mean(ConsRows$inccand), median(ConsRows$absdepdist), mean(LabRows$outsidecon))

xLD<-c(mean(LDRows$conservatives), mean(LDRows$liberal_democrats), mean(LDRows$partyfeeling), l, mean(LDRows$inccand), median(LDRows$absdepdist), mean(LabRows$outsidecon))


#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
		#numerator of pi 
		catL<-exp(xL%*%MIbetas[i,])
		catC<-exp(xC%*%MIbetas[i,])
		catLD<-exp(xLD%*%MIbetas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piL[i]<-catL/sum
		piC[i]<-catC/sum
		piLD[i]<-catLD/sum
	# Capturing mean Conservative vote share plus upper and lower 	bounds
	p.ests.cons[i,m]<-piC[i]
	}

}

#Without 50% line with new y-axis
pdf("ImputedDataLine.pdf")
plot(dists, apply(p.ests.cons,2,mean), ylim=c(.05, .65), col="white", main="Conservative vote share by distance with imputed data", xlab="Conservative candidate distance in km", ylab="Probability")
lines(dists, apply(p.ests.cons, 2, mean), col="blue", lwd=2)
lines(dists, apply(p.ests.cons, 2, quantile, .025, col="blue"), lty=2)
lines(dists, apply(p.ests.cons, 2, quantile, .975, col="blue"), lty=2)
#abline(h=.5, lty=2, col="red")
legend("topright", c("Mean", "95% confidence bounds"), lty=c(1, 2), lwd=c(3, 2), col=c("blue", "black"))
dev.off()
first.diffs<-p.ests.cons[,1]
min(first.diffs)


##-------------------------------------------------
#TERNARY PLOT

# Insert real mean distances for each party
		j<-mean(IMPUTE$drdistp)
		k<-120
		l<-mean(IMPUTE$drdistp)
	# Establish X values
xL<-c(mean(LabRows$conservatives), mean(LabRows$liberal_democrats), mean(LabRows$partyfeeling), j, mean(LabRows$inccand), median(LabRows$absdepdist), mean(LabRows$outsidecon))

xC<-c(mean(ConsRows$conservatives), mean(ConsRows$liberal_democrats), mean(ConsRows$partyfeeling), k, mean(ConsRows$inccand), median(ConsRows$absdepdist), mean(LabRows$outsidecon))

xLD<-c(mean(LDRows$conservatives), mean(LDRows$liberal_democrats), mean(LDRows$partyfeeling), l, mean(LDRows$inccand), median(LDRows$absdepdist), mean(LabRows$outsidecon))


#empty vectors for pi values to go in 
piL<-NULL
piC<-NULL
piLD<-NULL


#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
		#numerator of pi 
		catL<-exp(xL%*%MIbetas[i,])
		catC<-exp(xC%*%MIbetas[i,])
		catLD<-exp(xLD%*%MIbetas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piL[i]<-catL/sum
		piC[i]<-catC/sum
		piLD[i]<-catLD/sum
	}

mean(piL)
mean(piC)
mean(piLD)

mat<-cbind(piL, piC, piLD)
mean.point<-cbind(mean(piL), mean(piC), mean(piLD))
mat.2<-rbind(mat, mean.point)

# Coloring points by predicted probabilities of party vote
color<-matrix(nrow=nrow(mat), ncol=1)
for (i in 1:nrow(mat)){
	ifelse(max(mat[i,])==mat[i,1], color[i]<-"red", color[i]<-color[i])
	ifelse(max(mat[i,])==mat[i,2], color[i]<-"blue", color[i]<-color[i])
	ifelse(max(mat[i,])==mat[i,3], color[i]<-"yellow", color[i]<-color[i])
}
head(color)
color.2<-cbind(color, "black")


# Ternary graph!
pdf("NewTP.pdf")
ternaryplot(mat.2, scale=1, pch=16, cex=.3, col=color.2, dimnames=c("Labour", "Conservative", "Liberal Democrat"), dimnames_position="corner", main="New")
dev.off()

#------------------------
#simulate distances 
data<-read.csv("~/Desktop/PGreplication.csv")

head(data)
#class(data$aq16_)
#install.packages("survival") 
library(survival)

# Subset and clean data: have to change the variable that lists 0-strongly dislike instead of just using integers 
newdata<-subset(data, touse==1)
head(newdata)
newdata$aq16_=as.character(newdata$aq16_)
newdata$aq16_=ifelse(newdata$aq16_=="", NA, newdata$aq16_)
newdata$aq16_=ifelse(newdata$aq16_=="0-strongly dislike", 0, newdata$aq16_)
newdata$aq16_=ifelse(newdata$aq16_=="10-strongly like", 10, newdata$aq16_)
newdata$aq16_=as.numeric(newdata$aq16_)
newdata$aq16_=ifelse(newdata$aq16_=="0-strongly agree", 0, newdata$aq16_)

clogit1.4<-clogit(votechoice~conservatives+liberal_democrats+aq16_+drdistp+ inccand+absdepdist+strata(id)+cluster(aconsnam), method="approximate", data=newdata)
clogit1.4
mean4<-coef(clogit1.4)
vcv1.4<-vcov(clogit1.4)
library(MASS)
#there's an error here we don't understand 
set.seed(02138)
betas<-mvrnorm(1000, mean4, vcv1.4) 

LabRows<-subset(newdata, newdata$partyvote==1)

#Conservative matrix: 
ConsRows<-subset(newdata, newdata$partyvote==2)

#Liberal Democrat matrix: 
LDRows<-subset(newdata, newdata$partyvote==3)

dists<-seq(0, 120, by=5)

p.ests.conso<-matrix(data=NA, ncol=length(dists), nrow=1000)

for (m in 1:length(dists)){
# empty vectors for means and bounds

#empty vectors for pi values to go in 
piLo<-NULL
piCo<-NULL
piLDo<-NULL

# Insert real mean distances for each party
		j<-mean(newdata$drdistp, na.rm=T)
		k<-dists[m]
		l<-mean(newdata$drdistp, na.rm=T)
	# Establish X values
xL<-c(mean(LabRows$conservatives, na.rm=T), mean(LabRows$liberal_democrats, na.rm=T), mean(LabRows$aq16_, na.rm=T), j, mean(LabRows$inccand, na.rm=T), mean(LabRows$absdepdist, na.rm=T))

xC<-c(mean(ConsRows$conservatives, na.rm=T), mean(ConsRows$liberal_democrats, na.rm=T), mean(ConsRows$aq16_, na.rm=T), k, mean(ConsRows$inccand, na.rm=T), mean(ConsRows$absdepdist, na.rm=T))

xLD<-c(mean(LDRows$conservatives, na.rm=T), mean(LDRows$liberal_democrats, na.rm=T), mean(LDRows$aq16_, na.rm=T), l, mean(LDRows$inccand, na.rm=T), mean(LDRows$absdepdist, na.rm=T))


#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
		#numerator of pi 
		catL<-exp(xL%*%betas[i,])
		catC<-exp(xC%*%betas[i,])
		catLD<-exp(xLD%*%betas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piLo[i]<-catL/sum
		piCo[i]<-catC/sum
		piLDo[i]<-catLD/sum
		
	# Capturing mean Conservative vote share plus upper and lower 	bounds
	p.ests.conso[i,m]<-piCo[i]
}

}

plot(dists, apply(p.ests.conso,2,mean), ylim=c(.2, .65), col="white", main="Conservative vote share by distance", xlab="Conservative candidate distance in km", ylab="Probability")
lines(dists, apply(p.ests.conso, 2, mean), col="blue", lwd=2)
lines(dists, apply(p.ests.conso, 2, quantile, .025, col="blue"), lty=2)
lines(dists, apply(p.ests.conso, 2, quantile, .975, col="blue"), lty=2)
abline(h=.5, lty=2, col="red")
legend("topright", c("Mean", "95% confidence bounds", "50% probability"), lty=c(1, 2, 2), lwd=c(3, 2, 2), col=c("blue", "black", "red"))

# New y-axis
pdf("ConsNEW.pdf")
plot(dists, apply(p.ests.conso,2,mean), ylim=c(.05, .65), col="white", main="Conservative vote share by distance", xlab="Conservative candidate distance in km", ylab="Vote share")
lines(dists, apply(p.ests.conso, 2, mean), col="blue", lwd=2)
lines(dists, apply(p.ests.conso, 2, quantile, .025, col="blue"), col="blue", lty=2)
lines(dists, apply(p.ests.conso, 2, quantile, .975, col="blue"), col="blue", lty=2)
lines(dists, apply(p.ests.cons,2,mean), col="green4", lwd=2)
lines(dists, apply(p.ests.cons, 2, quantile, .025, col="green4"), col="green4", lty=2) 
lines(dists, apply(p.ests.cons, 2, quantile, .975, col="yellow"), col="green4", lty=2)  

legend("bottomright", c("Original Mean", "New Mean", "Old 95% confidence bounds", "New 95% confidence bounds"), lty=c(1, 1, 2, 2), lwd=c(3, 3, 2, 2), col=c("blue", "green4", "blue", "green4"))
dev.off()

###########---------------
#Labour confidence intervals
#randomly draw betas 

MIbetas<-matrix(nrow=1000, ncol=3)
MIbetas<-rbind(betasI1, betasI2, betasI3, betasI4, betasI5)

IMPUTE<-rbind(Idata1, Idata2, Idata3, Idata4, Idata5)
IMPUTE<-as.data.frame(IMPUTE)

####RESULTS for IMPUTATION
##Subset by Row
LabRows<-subset(IMPUTE, IMPUTE$partyvote==1)

#Conservative matrix: 
ConsRows<-subset(IMPUTE, IMPUTE$partyvote==2)

#Liberal Democrat matrix: 
LDRows<-subset(IMPUTE, IMPUTE$partyvote==3)


#simulate distances 
dists<-seq(0, 120, by=5)

p.ests.consLnew<-matrix(data=NA, ncol=length(dists), nrow=1000)

for (m in 1:length(dists)){

#empty vectors for pi values to go in 
piLL<-NULL
piCL<-NULL
piLDL<-NULL

# Insert real mean distances for each party
		j<-dists[m]
		k<-mean(IMPUTE$drdistp)
		l<-mean(IMPUTE$drdistp)
	# Establish X values
xL<-c(mean(LabRows$conservatives), mean(LabRows$liberal_democrats), mean(LabRows$partyfeeling), j, mean(LabRows$inccand), median(LabRows$absdepdist), mean(LabRows$outsidecon))

xC<-c(mean(ConsRows$conservatives), mean(ConsRows$liberal_democrats), mean(ConsRows$partyfeeling), k, mean(ConsRows$inccand), median(ConsRows$absdepdist), mean(LabRows$outsidecon))

xLD<-c(mean(LDRows$conservatives), mean(LDRows$liberal_democrats), mean(LDRows$partyfeeling), l, mean(LDRows$inccand), median(LDRows$absdepdist), mean(LabRows$outsidecon))


#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
		#numerator of pi 
		catL<-exp(xL%*%MIbetas[i,])
		catC<-exp(xC%*%MIbetas[i,])
		catLD<-exp(xLD%*%MIbetas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piLL[i]<-catL/sum
		piCL[i]<-catC/sum
		piLDL[i]<-catLD/sum
	# Capturing mean Conservative vote share plus upper and lower 	bounds
	p.ests.consLnew[i,m]<-piLL[i]
	}

}



#----------------------------
#original Labour vs. New Labour

set.seed(02138)
betas<-mvrnorm(1000, mean4, vcv1.4) 

LabRows<-subset(newdata, newdata$partyvote==1)

#Conservative matrix: 
ConsRows<-subset(newdata, newdata$partyvote==2)

#Liberal Democrat matrix: 
LDRows<-subset(newdata, newdata$partyvote==3)

dists<-seq(0, 120, by=5)

p.ests.consLoff<-matrix(data=NA, ncol=length(dists), nrow=1000)

for (m in 1:length(dists)){
# empty vectors for means and bounds

#empty vectors for pi values to go in 
piLLod<-NULL
piCLod<-NULL
piLDLod<-NULL

# Insert real mean distances for each party
		j<-dists[m]
		k<-mean(newdata$drdistp, na.rm=T)
		l<-mean(newdata$drdistp, na.rm=T)
	# Establish X values
xL<-c(mean(LabRows$conservatives, na.rm=T), mean(LabRows$liberal_democrats, na.rm=T), mean(LabRows$aq16_, na.rm=T), j, mean(LabRows$inccand, na.rm=T), mean(LabRows$absdepdist, na.rm=T))

xC<-c(mean(ConsRows$conservatives, na.rm=T), mean(ConsRows$liberal_democrats, na.rm=T), mean(ConsRows$aq16_, na.rm=T), k, mean(ConsRows$inccand, na.rm=T), mean(ConsRows$absdepdist, na.rm=T))

xLD<-c(mean(LDRows$conservatives, na.rm=T), mean(LDRows$liberal_democrats, na.rm=T), mean(LDRows$aq16_, na.rm=T), l, mean(LDRows$inccand, na.rm=T), mean(LDRows$absdepdist, na.rm=T))


#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
		#numerator of pi 
		catL<-exp(xL%*%betas[i,])
		catC<-exp(xC%*%betas[i,])
		catLD<-exp(xLD%*%betas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piLLod[i]<-catL/sum
		piCLod[i]<-catC/sum
		piLDLod[i]<-catLD/sum
		
	# Capturing mean Conservative vote share plus upper and lower 	bounds
	p.ests.consLoff[i,m]<-piLLod[i]
}

}

#Without 50% line with new y-axis
pdf("ImputedDataLineL.pdf")
plot(dists, apply(p.ests.consL,2,mean), ylim=c(.05, .65), col="white", main="Labour vote share by distance", xlab="Labour candidate distance in km", ylab="Vote share")
lines(dists, apply(p.ests.consLnew, 2, mean), col="purple4", lwd=2)
lines(dists, apply(p.ests.consLnew, 2, quantile, .025, col="blue"), col="purple4", lty=2)
lines(dists, apply(p.ests.consLnew, 2, quantile, .975, col="blue"), col="purple4", lty=2)
lines(dists, apply(p.ests.consLoff, 2, mean), col="red", lwd=2)
lines(dists, apply(p.ests.consLoff, 2, quantile, .025, col="blue"), col="red", lty=2)
lines(dists, apply(p.ests.consLoff, 2, quantile, .975, col="blue"), col="red", lty=2)
#abline(h=.5, lty=2, col="red")
legend("topright", c("Original Mean", "New Mean", "Old 95% confidence bounds", "New 95% confidence bounds"), lty=c(1, 1, 2, 2), lwd=c(3, 3, 2, 2), col=c("red", "purple4", "red", "purple4"))
dev.off()




###########---------------
#LD New
#randomly draw betas 

MIbetas<-matrix(nrow=1000, ncol=3)
MIbetas<-rbind(betasI1, betasI2, betasI3, betasI4, betasI5)

IMPUTE<-rbind(Idata1, Idata2, Idata3, Idata4, Idata5)
IMPUTE<-as.data.frame(IMPUTE)

####RESULTS for IMPUTATION
##Subset by Row
LabRows<-subset(IMPUTE, IMPUTE$partyvote==1)

#Conservative matrix: 
ConsRows<-subset(IMPUTE, IMPUTE$partyvote==2)

#Liberal Democrat matrix: 
LDRows<-subset(IMPUTE, IMPUTE$partyvote==3)


#simulate distances 
dists<-seq(0, 120, by=5)

p.ests.consLoo<-matrix(data=NA, ncol=length(dists), nrow=1000)

for (m in 1:length(dists)){

#empty vectors for pi values to go in 
piLLoo<-NULL
piCLoo<-NULL
piLDLoo<-NULL

# Insert real mean distances for each party
		j<-mean(IMPUTE$drdistp)
		k<-mean(IMPUTE$drdistp)
		l<-dists[m]
	# Establish X values
xL<-c(mean(LabRows$conservatives), mean(LabRows$liberal_democrats), mean(LabRows$partyfeeling), j, mean(LabRows$inccand), median(LabRows$absdepdist), mean(LabRows$outsidecon))

xC<-c(mean(ConsRows$conservatives), mean(ConsRows$liberal_democrats), mean(ConsRows$partyfeeling), k, mean(ConsRows$inccand), median(ConsRows$absdepdist), mean(LabRows$outsidecon))

xLD<-c(mean(LDRows$conservatives), mean(LDRows$liberal_democrats), mean(LDRows$partyfeeling), l, mean(LDRows$inccand), median(LDRows$absdepdist), mean(LabRows$outsidecon))


#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
		#numerator of pi 
		catL<-exp(xL%*%MIbetas[i,])
		catC<-exp(xC%*%MIbetas[i,])
		catLD<-exp(xLD%*%MIbetas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piLLoo[i]<-catL/sum
		piCLoo[i]<-catC/sum
		piLDoo[i]<-catLD/sum
	# Capturing mean Conservative vote share plus upper and lower 	bounds
	p.ests.consLoo[i,m]<-piLDoo[i]
	}

}



#----------------------------
#original LD

set.seed(02138)
betas<-mvrnorm(1000, mean4, vcv1.4) 

LabRows<-subset(newdata, newdata$partyvote==1)

#Conservative matrix: 
ConsRows<-subset(newdata, newdata$partyvote==2)

#Liberal Democrat matrix: 
LDRows<-subset(newdata, newdata$partyvote==3)

dists<-seq(0, 120, by=5)

p.ests.consLod<-matrix(data=NA, ncol=length(dists), nrow=1000)

for (m in 1:length(dists)){
# empty vectors for means and bounds

#empty vectors for pi values to go in 
piLod<-NULL
piCLod<-NULL
piLDod<-NULL

# Insert real mean distances for each party
		j<-mean(newdata$drdistp, na.rm=T)
		k<-mean(newdata$drdistp, na.rm=T)
		l<-dists[m]
	# Establish X values
xL<-c(mean(LabRows$conservatives, na.rm=T), mean(LabRows$liberal_democrats, na.rm=T), mean(LabRows$aq16_, na.rm=T), j, mean(LabRows$inccand, na.rm=T), mean(LabRows$absdepdist, na.rm=T))

xC<-c(mean(ConsRows$conservatives, na.rm=T), mean(ConsRows$liberal_democrats, na.rm=T), mean(ConsRows$aq16_, na.rm=T), k, mean(ConsRows$inccand, na.rm=T), mean(ConsRows$absdepdist, na.rm=T))

xLD<-c(mean(LDRows$conservatives, na.rm=T), mean(LDRows$liberal_democrats, na.rm=T), mean(LDRows$aq16_, na.rm=T), l, mean(LDRows$inccand, na.rm=T), mean(LDRows$absdepdist, na.rm=T))


#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
		#numerator of pi 
		catL<-exp(xL%*%betas[i,])
		catC<-exp(xC%*%betas[i,])
		catLD<-exp(xLD%*%betas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piLod[i]<-catL/sum
		piCLod[i]<-catC/sum
		piLDod[i]<-catLD/sum
		
	# Capturing mean Conservative vote share plus upper and lower 	bounds
	p.ests.consLod[i,m]<-piLDod[i]
}

}

#Without 50% line with new y-axis
pdf("ImputedDataLineLD.pdf")
plot(dists, apply(p.ests.consL,2,mean), ylim=c(.05, .65), col="white", main="Liberal Democratic vote share by distance", xlab="Liberal Democratic candidate distance in km", ylab="Vote share")
lines(dists, apply(p.ests.consLoo, 2, mean), col="hotpink", lwd=2)
lines(dists, apply(p.ests.consLoo, 2, quantile, .025, col="blue"), col="hotpink", lty=2)
lines(dists, apply(p.ests.consLoo, 2, quantile, .975, col="blue"), col="hotpink", lty=2)
lines(dists, apply(p.ests.consLod, 2, mean), col="goldenrod", lwd=2)
lines(dists, apply(p.ests.consLod, 2, quantile, .025, col="blue"), col="goldenrod", lty=2)
lines(dists, apply(p.ests.consLod, 2, quantile, .975, col="blue"), col="goldenrod", lty=2)
#abline(h=.5, lty=2, col="red")
legend("topright", c("Original Mean", "New Mean", "Old 95% confidence bounds", "New 95% confidence bounds"), lty=c(1, 1, 2, 2), lwd=c(3, 3, 2, 2), col=c("goldenrod", "hotpink", "goldenrod", "hotpink"))
dev.off()
